function [outputresultmap outputresulterrormap] = ratheartregister(fixedstudydir, movingstudydir)
% ratheartregister: run registration between multiple days
%
% Usage:  outputresultmap = ratheartregister(fixedstudydir, movingstudydir)
%
% Where: 
%   fixedstudydir, movingstudydir are the T2 directories of day 1, and
% a next day, respectively.
%   outputresultmap = the name of the registered (moving) R2map dcm.
%

% remove slashes at the end: fixed
while ( fixedstudydir(end) == '/' || fixedstudydir(end) == '\' )
    fixedstudydir = fixedstudydir(1:(end-1));
end
% output dir
fixeddir = [ fixedstudydir 'postproc\' ];
if ~exist(fixeddir, 'dir')
  error( 'ratheartregister:invalidoutputdir', 'Invalid postproc dir. Run first ratheartT2fit!' );  
end

% remove slashes at the end: moving
while ( movingstudydir(end) == '/' || movingstudydir(end) == '\' )
    movingstudydir = movingstudydir(1:(end-1));
end
% output dir
movingdir = [ movingstudydir 'postproc\' ];
if ~exist(movingdir, 'dir')
  error( 'ratheartregister:invalidoutputdir', 'Invalid postproc dir. Run first ratheartT2fit!' );  
end

% find elastix.exe
currentdir= fileparts( which('dGEmricIRregister.m') );
cd(currentdir)
cd ..
elastixdir = [pwd filesep 'elastix'];

% define all input arguments for elastix
fixedim = fullfile(fixeddir, 'FirstIm.dcm');
movingim = fullfile(movingdir, 'FirstIm.dcm');
fixedmask = ratheartmakeregmask(fixedstudydir);
movingmask = ratheartmakeregmask(movingstudydir);
fixedpoints = fullfile(fixeddir, 'points.txt');
movingpoints = fullfile(movingdir, 'points.txt');
par1 = fullfile(elastixdir, 'parameters_RigidPointBased.txt');
par2 = fullfile(elastixdir, 'parameters_RigidIntensityBased.txt');
outdir = fullfile(movingdir, 'reg');
finaltp = fullfile(outdir, 'TransformParameters.1.txt');
resultmap = fullfile(outdir, 'result.dcm');
resultTE = fullfile(outdir, 'result.1.dcm');

outputresultTE = fullfile(movingdir, 'registeredFirstIm.dcm');

%define R1/R2 map
cd(movingdir)
d = dir;
for j=3:size(d,1)
    x=d(j).name;
    if size(char(regexp(x,'R1map', 'match')),2)>=1 %read in the fixed longest IR image
        movingmap = fullfile(movingdir, 'R1map.dcm');
        movingerrormap = fullfile(movingdir, 'R1errormap.dcm');
        outputresultmap = fullfile(movingdir, 'registeredR1map.dcm');
        outputresulterrormap = fullfile(movingdir, 'registeredR1errormap.dcm');
        movingmapT = fullfile(movingdir, 'T1map.dcm');
        movingerrormapT = fullfile(movingdir, 'T1errormap.dcm');
        outputresultmapT = fullfile(movingdir, 'registeredT1map.dcm');
        outputresulterrormapT = fullfile(movingdir, 'registeredT1errormap.dcm');
    elseif size(char(regexp(x,'R2map', 'match')),2)>=1 %read in the fixed longest IR image
        movingmap = fullfile(movingdir, 'R2map.dcm');
        movingerrormap = fullfile(movingdir, 'R2errormap.dcm');
        outputresultmap = fullfile(movingdir, 'registeredR2map.dcm');
        outputresulterrormap = fullfile(movingdir, 'registeredR2errormap.dcm');
        movingmapT = fullfile(movingdir, 'T2map.dcm');
        movingerrormapT = fullfile(movingdir, 'T2errormap.dcm');
        outputresultmapT = fullfile(movingdir, 'registeredT2map.dcm');
        outputresulterrormapT = fullfile(movingdir, 'registeredT2errormap.dcm');
    end   
end

if ~exist(outdir, 'dir')
  disp( ['making ' outdir]);     
  [success, message, messageid] = mkdir(outdir);
  if ~success    
     error( 'ratheartregister:invalidoutputdir', message );     
  end 
end

checkfiles = {fixedim, movingim, fixedmask, movingmask, fixedpoints, movingpoints, par1, par2, movingmap};
for cf = 1:numel(checkfiles)
  if ~exist(checkfiles{cf})
      disp( [checkfiles{cf} ' is missing!' ] );
      error( 'ratheartregister:missinginputfile', 'missing input file!' );
  end
end

% call elastix 
cmdline1 = ['"' fullfile(elastixdir, 'elastix') ...
  '" -f "' fixedim '" -m "' movingim ...
  '" -fp "' fixedpoints '" -mp "' movingpoints ...
  '" -fMask "' fixedmask '" -mMask "' movingmask ...
  '" -p "' par1 '" -p "' par2 '" -out "' outdir '"'];

[status1,result1] = system(cmdline1);

% check if elastix was successful
checkfiles = {finaltp};
for cf = 1:numel(checkfiles)
  if ~exist(checkfiles{cf})
      disp( [checkfiles{cf} ' is missing!' ] );
      disp( result1 );
      error( 'ratheartregister:missingoutputfile', 'missing output file from elastix!' );      
  end
end

% copy the result.1.mhd to outputresultTE
copyfile( resultTE, outputresultTE );

% call transformix, to transform the R2map
cmdline2 = ['"' fullfile(elastixdir, 'transformix') '" -in "' movingmap '" -out "' outdir '" -tp "' finaltp '"'];
[status2,result2] = system(cmdline2);

% copy the result to outputresultmap
copyfile( resultmap, outputresultmap );

% call transformix, to transform the R2errormap
cmdline2 = ['"' fullfile(elastixdir, 'transformix') '" -in "' movingerrormap '" -out "' outdir '" -tp "' finaltp '"'];
[status2,result2] = system(cmdline2);

% copy the result to outputresultmap
copyfile( resultmap, outputresulterrormap );

% call transformix, to transform the T2map
cmdline2 = ['"' fullfile(elastixdir, 'transformix') '" -in "' movingmapT '" -out "' outdir '" -tp "' finaltp '"'];
[status2,result2] = system(cmdline2);

% copy the result to outputresultmap
copyfile( resultmap, outputresultmapT );

% call transformix, to transform the T2errormap
cmdline2 = ['"' fullfile(elastixdir, 'transformix') '" -in "' movingerrormapT '" -out "' outdir '" -tp "' finaltp '"'];
[status2,result2] = system(cmdline2);

% copy the result to outputresultmap
copyfile( resultmap, outputresulterrormapT );

